///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/vector/vector.h
///
/// \brief   header file for the formic::vector class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef FORMIC_VECTOR_VECTOR_HEADER
#define FORMIC_VECTOR_VECTOR_HEADER

#include <iostream>
#include <cassert>
#include <algorithm>
//#include <list>
#include <map>
#include <stack>
#include <vector>

#include <boost/shared_ptr.hpp>
#include <boost/shared_array.hpp>

#include <formic/exception.h>

namespace formic {

  template<class ScalarType, class LengthType> class VectorLock {

    private:

      typedef ScalarType S;
      typedef LengthType L;

    private:

      L _length;

      boost::shared_array<S> _array;

      //static std::list<L>                       available_lengths;

      //static std::list<boost::shared_array<S> > available_arrays;

      static std::map<L, std::stack<boost::shared_array<S> > > available_arrays;

    private:

      VectorLock & operator=(const VectorLock &);

      VectorLock(const VectorLock &);

      static void pthread_lock() {}

      static void pthread_unlock() {}

    public:

      L length() { return _length; }

      S * cptr() { return _array.get(); }

      VectorLock(const L n);

      ~VectorLock();

      static void garbage_collect();

  };

  void vector_garbage_collect();


  template<class ScalarType, class LengthType> class Matrix;
  template<class ScalarType, class LengthType> class ConstMatrix;


  template<class ScalarType, class LengthType> class Vector {

    private:

      typedef ScalarType S;
      typedef LengthType L;

    private:

      L _n;

      S * _p;

      boost::shared_ptr<formic::VectorLock<S,L> > _lock;

    private:

      void initialize_with_new_lock(const L n) {
        if (n < 0)
          throw formic::Exception("cannot initialize a vector with a negative length");
        _n = n;
        _p = 0;
        _lock.reset();
        if ( n > 0 ) {
          boost::shared_ptr<formic::VectorLock<S,L> > new_lock( new formic::VectorLock<S,L>(n) );
          _lock = new_lock;
          _p = _lock->cptr();
        }
      }

    public:

      // constructor without initialization of elements
      explicit Vector(const L n = 0) {
        this->initialize_with_new_lock(n);
      }

      // constructor with initialization of elements
      explicit Vector(const L n, const S & v) {
        this->initialize_with_new_lock(n);
        if (_n > 0)
          std::fill(_p, _p + _n, v);
      }

      // copy constructer (shallow copy)
      Vector(const Vector & other) {
        _n = other._n;
        _p = other._p;
        _lock = other._lock;
      }

      // assignment operator (shallow copy)
      Vector & operator=(const Vector & other) {
        //std::cout << boost::format("assigning size %i") % other._n << std::endl;
        if (this != &other) {
          _n = other._n;
          _p = other._p;
          _lock = other._lock;
        }
        //std::cout << boost::format("assignment complete")<< std::endl;
        return *this;
      }

      // deep copy to c-style array
      const Vector & operator>>=(S * const x) const {
        std::copy(_p, _p + _n, x);
        return *this;
      }

      // deep copy from c-style array
      Vector & operator<<=(const S * const x) {
        std::copy(x, x + _n, _p);
        return *this;
      }

      // deep copy to std::vector
      const Vector & operator>>=(std::vector<S> & x) const {
        if ( this->size() != x.size() )
          x.resize(this->size());
        if ( this->size() > 0 )
          *this >>= &x.at(0);
        return *this;
      }

      // deep copy from std::vector
      Vector & operator<<=(const std::vector<S> & x) {
        if ( this->size() != x.size() )
          this->reset(x.size());
        if ( x.size() > 0 )
          *this <<= &x.at(0);
        return *this;
      }

      // deep copy assignment
      Vector & operator<<=(const Vector & other) {
        if (this != &other) {
          if ( this->size() != other.size() )
            this->initialize_with_new_lock(other._n);
          *this <<= other._p;
        }
        return *this;
      }

      /// \brief  clears the vector, leaving it as a 0 length vector
      void clear() { *this = formic::Vector<S,L>(); }

      /// \brief  clears the vector and re-allocates it to be the specified size
      Vector & reset(const int n) {
        this->clear();
        return *this = formic::Vector<S,L>(n);
      }

      /// \brief  clears the vector and re-allocates it to be the specified size, initializing its elements to v
      Vector & reset(const int n, const S & v) {
        this->clear();
        return *this = formic::Vector<S,L>(n,v);
      }

      ///// \brief  clears the vector and re-allocates it to be the specified size, initializing its elements with the array pointed to by x
      //Vector & reset(const int n, const S * const x) {
      //  return this->reset(n) <<= x;
      //}

      const S & operator()(L i) const {
        #ifndef NDEBUG
        if ( i < 0 || i >= _n ) throw formic::Exception("out-of-range index i = %i in formic::Vector::operator() for vector of length %i") % i % _n;
        #endif
        return _p[i];
      }

      S & operator()(L i) {
        #ifndef NDEBUG
        if ( i < 0 || i >= _n ) throw formic::Exception("out-of-range index i = %i in formic::Vector::operator() for vector of length %i") % i % _n;
        #endif
        return _p[i];
      }

      const S & operator[](L i) const {
        #ifndef NDEBUG
        if ( i < 0 || i >= _n ) throw formic::Exception("out-of-range index i = %i in formic::Vector::operator[] for vector of length %i") % i % _n;
        #endif
        return _p[i];
      }

      S & operator[](L i) {
        #ifndef NDEBUG
        if ( i < 0 || i >= _n ) throw formic::Exception("out-of-range index i = %i in formic::Vector::operator[] for vector of length %i") % i % _n;
        #endif
        return _p[i];
      }

      const S & at(L i) const {
        if ( i < 0 || i >= _n ) throw formic::Exception("out-of-range index i = %i in formic::Vector::at for vector of length %i") % i % _n;
        return _p[i];
      }

      S & at(L i) {
        if ( i < 0 || i >= _n ) throw formic::Exception("out-of-range index i = %i in formic::Vector::at for vector of length %i") % i % _n;
        return _p[i];
      }

      S *       begin()       { return _p; }
      const S * begin() const { return _p; }

      S *       end()         { return _p + this->size(); }
      const S * end()   const { return _p + this->size(); }

      formic::Vector<S,L> clone() const {
        formic::Vector<S,L> retval;
        retval <<= *this;
        return retval;
      }

      Vector & operator=(const S a) {
        //for (L i = 0; i < _n; i++)
        //  _p[i] = a;
        if (_n > 0)
          std::fill(_p, _p + _n, a);
        return *this;
      }

      formic::Vector<S,L> & operator+=(const S a) {
        for (L i = 0; i < _n; i++)
          _p[i] += a;
        return *this;
      }

      formic::Vector<S,L> & operator-=(const S a) {
        for (L i = 0; i < _n; i++)
          _p[i] -= a;
        return *this;
      }

      formic::Vector<S,L> & operator*=(const S a) {
        for (L i = 0; i < _n; i++)
          _p[i] *= a;
        return *this;
      }

      formic::Vector<S,L> & operator/=(const S a) {
        for (L i = 0; i < _n; i++)
          _p[i] /= a;
        return *this;
      }

      formic::Vector<S,L> & operator+=(const formic::Vector<S,L> & other) {
        if (_n != other._n)
          throw formic::Exception("cannot add vectors of different lengths");
        for (L i = 0; i < _n; i++)
          _p[i] += other._p[i];
        return *this;
      }

      formic::Vector<S,L> & operator-=(const formic::Vector<S,L> & other) {
        if (_n != other._n)
          throw formic::Exception("cannot subtract vectors of different lengths");
        for (L i = 0; i < _n; i++)
          _p[i] -= other._p[i];
        return *this;
      }

      formic::Vector<S,L> & operator*=(const formic::Vector<S,L> & other) {
        if (_n != other._n)
          throw formic::Exception("cannot multiply vectors of different lengths");
        for (L i = 0; i < _n; i++)
          _p[i] *= other._p[i];
        return *this;
      }

      formic::Vector<S,L> & operator/=(const formic::Vector<S,L> & other) {
        if (_n != other._n)
          throw formic::Exception("cannot divide vectors of different lengths");
        for (L i = 0; i < _n; i++)
          _p[i] /= other._p[i];
        return *this;
      }

      formic::Vector<S,L> & axpy(const S a, const formic::Vector<S,L> & other) {
        if (_n != other._n)
          throw formic::Exception("cannot axpy vectors of different lengths");
        for (L i = 0; i < _n; i++)
          _p[i] += a * other._p[i];
        return *this;
      }

      formic::Vector<S,L> & axpy(const S a, const S * const x) {
        for (L i = 0; i < _n; i++)
          _p[i] += a * x[i];
        return *this;
      }

      formic::Vector<S,L> operator-() {
        formic::Vector<S,L> retval(_n);
        for (L i = 0; i < _n; i++)
          retval._p[i] = -_p[i];
        return retval;
      }

      formic::Vector<S,L> operator+() {
        return this->clone();
      }

      L size() const { return _n; }

      formic::Matrix<S,L> t() const;
      formic::Matrix<S,L> m() const;

  };

}

template<class S, class L>
formic::Vector<S,L> operator+(const formic::Vector<S,L> & v1, const formic::Vector<S,L> & v2) {
  if ( v1.size() != v2.size() )
    throw formic::Exception("cannot add vectors of different lengths");
  //return ( formic::Vector<S,L>(v1) += v2 );
  formic::Vector<S,L> retval = v1.clone();
  retval += v2;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator-(const formic::Vector<S,L> & v1, const formic::Vector<S,L> & v2) {
  if ( v1.size() != v2.size() )
    throw formic::Exception("cannot subtract vectors of different lengths");
  formic::Vector<S,L> retval = v1.clone();
  retval -= v2;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator*(const formic::Vector<S,L> & v1, const formic::Vector<S,L> & v2) {
  if ( v1.size() != v2.size() )
    throw formic::Exception("cannot multiply vectors of different lengths");
  formic::Vector<S,L> retval = v1.clone();
  retval *= v2;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator/(const formic::Vector<S,L> & v1, const formic::Vector<S,L> & v2) {
  if ( v1.size() != v2.size() )
    throw formic::Exception("cannot divide vectors of different lengths");
  formic::Vector<S,L> retval = v1.clone();
  retval /= v2;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator+(const formic::Vector<S,L> & v1, const S a) {
  formic::Vector<S,L> retval = v1.clone();
  retval += a;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator-(const formic::Vector<S,L> & v1, const S a) {
  formic::Vector<S,L> retval = v1.clone();
  retval -= a;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator*(const formic::Vector<S,L> & v1, const S a) {
  formic::Vector<S,L> retval = v1.clone();
  retval *= a;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator/(const formic::Vector<S,L> & v1, const S a) {
  formic::Vector<S,L> retval = v1.clone();
  retval /= a;
  return retval;
}




template<class S, class L>
formic::Vector<S,L> operator+(const S a, const formic::Vector<S,L> & v1) {
  formic::Vector<S,L> retval = v1.clone();
  retval += a;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator-(const S a, const formic::Vector<S,L> & v1) {
  formic::Vector<S,L> retval(v1.size(), a);
  retval -= v1;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator*(const S a, const formic::Vector<S,L> & v1) {
  formic::Vector<S,L> retval = v1.clone();
  retval *= a;
  return retval;
}

template<class S, class L>
formic::Vector<S,L> operator/(const S a, const formic::Vector<S,L> & v1) {
  formic::Vector<S,L> retval(v1.size(), a);
  retval /= v1;
  return retval;
}

#endif
